Consignes

Complétez ce document en remplissant les chunks vides pour écrire le code qui vous a permis de répondre à la question. Les réponses attendant un résultat chiffré ou une explication devront être insérés entre le balises html code. Par exemple pour répondre à la question suivante :

La bioinfo c'est : <code>MERVEILLEUX</code>.

N’hésitez pas à commenter votre code, enrichier le rapport en y insérant des résultats ou des graphiques/images pour expliquer votre démarche. N’oubliez pas les bonnes pratiques pour une recherche reproductible ! Nous souhaitons à minima que l’analyse soit reproductible sur le cluster de l’IFB.

Introduction

Vous allez travailler sur des données de reséquençage d’un génome bactérien : Bacillus subtilis. Les données sont issues de cet article :

Analyses

Organisation de votre espace de travail

## Repertoire de travail
pwd
#/shared/projects/dubii2021/agodmer

## Enregistrement de l'arborescence dans le fichier Evaluation_M4_M5/supplemental_data/Organisation_espace_de_travail.txt
tree > supplemental_data/Organisation_espace_de_travail.txt

## Commande tree
tree

#[agodmer@clust-slurm-client Evaluation_M4_M5]$ tree
#.
#|-- DATA_analysis
#|   |-- CLEANING
#|   |-- FASTQ
#|   |-- MAPPING
#|   |-- QC
#|   `-- REFERENCE_GENOME
#|-- EvaluationM4M5-main-results
#|   |-- Evaluation.Rmd
#|   |-- Evaluation.html
#|   |-- EvaluationM4M5.Rproj
#|   |-- README.md
#|   |-- css
#|   |   `-- style.css
#|   |-- images
#|   |   |-- inrae.png
#|   |   `-- migale-orange.png
#|   `-- resources
#|       |-- biblio.bib
#|       |-- biomed-central.csl
#|       `-- footer.html
#|-- LICENSE
#|-- README.md
#`-- supplemental_data
#    |-- Organisation_espace_de_travail.txt
#    `-- README_supp_data.md

 
## Il n'y a pas de modification à faire pour l'instant

Téléchargement des données brutes

Récupérez les fichiers FASTQ issus du run SRR10390685 grâce à l’outil sra-tools [1]

## Répertoire de travail
pwd
#/shared/projects/dubii2021/agodmer/Evaluation_M4_M5/EvaluationM4M5-main-results

## Charger le module SRA tools avec la dernière version 
module avail sra-tools
module load sra-tools/2.10.3

## Utilisation de la commande fasterq-dump pour télécharger les fichiers 

## Visualisation de la version
fasterq-dump --version
#"fasterq-dump" version 2.10.3

## Ecriture de la version de fasterq-dump sur le fichier Version_tools.txt dans supplemental_data
fasterq-dump --version > ../supplemental_data/Version_tools.txt

## Changement du répertoire de travail
cd ..

## Reservation des ressources pour le cluster
salloc --cpus-per-task=6 --mem=5G

## Télécharger les fichiers FASTQ
srun fasterq-dump --split-files -p SRR10390685 --outdir DATA_analysis/FASTQ/

## Lister les fichiers et regarder leur taille
ls -sh DATA_analysis/FASTQ/
#total 5.0G
#2.5G SRR10390685_1.fastq  2.5G SRR10390685_2.fastq

## Compression des fichier FASTQ avec mode verbeux acitvé
srun gzip --verbose DATA_analysis/FASTQ/*.fastq

# Liste des fichier et taille
ls -sh DATA_analysis/FASTQ/
#total 1.3G
#617M SRR10390685_1.fastq.gz  627M SRR10390685_2.fastq.gz

## La compression a bien fonctionnée !

Combien de reads sont présents dans les fichiers R1 et R2 ?

## Répertoire de travail
pwd
#/shared/projects/dubii2021/agodmer/Evaluation_M4_M5/DATA_analysis

## Chargement du module seqkit
module avail seqkit
module load seqkit/0.14.0
## Ecriture de la version de seqkit sur le fichier Version_tools.txt dans supplemental_data
seqkit version >> ../supplemental_data/Version_tools.txt

## Comptage des read présents et écriture des statistiques dans le fichier Raw_stats_fastq.txt dans Results/supplemental_data
srun seqkit stats --threads 1 FASTQ/*.fastq.gz > ../supplemental_data/Raw_stats_fastq.txt

## Visualisation des résultats
cat ../supplemental_data/Raw_stats_fastq.txt

#file                          format  type   num_seqs        sum_len  min_len  avg_len  max_len
#FASTQ/SRR10390685_1.fastq.gz  FASTQ   DNA   7,066,055  1,056,334,498       35    149.5      151
#FASTQ/SRR10390685_2.fastq.gz  FASTQ   DNA   7,066,055  1,062,807,718      130    150.4      151

Les fichiers FASTQ contiennent 7066055 reads.

Téléchargez le génome de référence de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse

## Repertoire de travail
pwd
#/shared/projects/dubii2021/agodmer/Evaluation_M4_M5/DATA_analysis

## Téléchargement du génome de référence
srun wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gz

## Visualisation du fichier
ls -sh
#total 68K
#68K GCF_000009045.1_ASM904v1_genomic.fna.gz

Quelle est la taille de ce génome ?


## Taille du génome de référence
## Enregistrement du résultats dans Stats_genome_ref.txt
srun zcat GCF_000009045.1_ASM904v1_genomic.fna.gz | grep -v "^>" | tr --delete "\n" | wc -c > ../../supplemental_data/Stats_genome_ref.txt

## Visualisation du fichier
cat ../../supplemental_data/Stats_genome_ref.txt
#4215606

La taille de ce génome est de 4215606 paires de bases.

Téléchargez l’annotation de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse

## Répertoire de travail
pwd
#/shared/projects/dubii2021/agodmer/Evaluation_M4_M5/DATA_analysis/REFERENCE_GENOME

## Téléchargement de l'annotation du génome
srun wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.gff.gz

## Visualisation des fichiers
ls -sh

#total 1.3M
#1.2M GCF_000009045.1_ASM904v1_genomic.fna.gz   69K GCF_000009045.1_ASM904v1_genomic.gff.gz

Combien de gènes sont connus pour ce génome ?

## Répertoire de travail
pwd
#/shared/projects/dubii2021/agodmer/Evaluation_M4_M5/DATA_analysis/REFERENCE_GENOME

## Enregistrement du résultats dans Stats_genome_ref_nb_genes.txt
zless GCF_000009045.1_ASM904v1_genomic.gff.gz | awk '($3 == "gene")' | wc -l > ../../supplemental_data/Stats_genome_ref_nb_genes.txt

## Visualisation du fichier
cat ../../supplemental_data/Stats_genome_ref_nb_genes.txt
#4448

4448 gènes sont recensés dans le fichier d’annotation.

Contrôle qualité

Lancez l’outil fastqc [2] dédié à l’analyse de la qualité des bases issues d’un séquençage haut-débit

## Répertoire de travail
pwd
#/shared/projects/dubii2021/agodmer/Evaluation_M4_M5/DATA_analysis

module avail fastqc
module load fasqc/0.11.9
## Ecriture de la version de fastqc sur le fichier Version_tools.txt dans supplemental_data
fastqc --version >> ../supplemental_data/Version_tools.txt

## Réservation des resources pour le cluster
salloc --cpus-per-task=8 --mem=10G

## Lancement de la commande fastqc
srun fastqc FASTQ/SRR10390685_1.fastq.gz -o QC/ -t 8
srun fastqc FASTQ/SRR10390685_2.fastq.gz -o QC/ -t 8

## Copie des rapport htlm vers supplemental data
scp QC/*html ../supplemental_data/

## Création d'un rapport multiQC

## Chargement du module multiQC
module avail multiqc
module load multiqc/1.9

## Ecriture de la version de multiqc avec nano sur le fichier Version_tools.txt
srun multiqc -d . -o .

## Copie du rapport dans supplmental_data
scp multiqc_report.html ../supplemental_data/

La qualité des bases vous paraît-elle satisfaisante ? Pourquoi ?

  • [X] Oui
  • [ ] Non

car la qualité (estimée par le score Phred) de chaque base pour tous les reads de votre jeu de données est
en moyenne et en médiane supérieur à 30 (cette valeur correspond à la probabilité de 1 pour 1000 q’une base ait été
identifiée incorectement) comme le montre le graphique Per base sequence quality qui représente représente
la qualité (score Phred, en ordonnée) de chaque base (en abscisse) pour tous les reads de votre jeu de données ;
Le code couleur vous indique les scores de très bonne qualité en vert pour la majorité de nos reads.
La qualité des reads sur la plupart des plates-formes se dégradera au fur et à mesure de la progression du séquençage,
il est donc courant de voir la qualité des reads entrer dans la zone orange vers la fin d’une lecture.
On peut observer un warning sur Per base sequence content qui trace le pourcentage de chacun des quatre nucléotides
à chaque position sur toutes les lectures dans le fichier de séquence d’entrée.
Cependant La proportion de chacune des quatre bases reste relativement constante sur la longueur de la lecture avec AT et GC
et les lignes du graphiques sont parallèles les unes aux autres.
La sur-représentation des bases AT est concordante avec le poucentage de GC de l’espèce étudiée (Bacillus subtilis) qui n’est pas égal
à 50% mais à 43% sur le génome de référence Bacillus subtilis subsp. subtilis str. 168 NCBI.

Est-ce que les reads déposés ont subi une étape de nettoyage avant d’être déposés ? Pourquoi ?

  • [X] Oui
  • [ ] Non

car Les reads n’ont pas la même longueur.
De plus ,on remarque que les adaptateurs ont été enlevés sur le graphique Adapter Content.

Quelle est la profondeur de séquençage (calculée par rapport à la taille du génome de référence) ?

## Répertoire de travail
pwd
#/shared/projects/dubii2021/agodmer/Evaluation_M4_M5

## Formule :  N*L/G ; avec N = Nombre de lectures, L = Longueur des lectures et G = Taille du génome

## taille du génome de référence (G)
zcat DATA_analysis/REFERENCE_GENOME/GCF_000009045.1_ASM904v1_genomic.fna.gz |grep -v "^>"|wc -c
#4268302

## Récupération des données : Nombre de lecture (N), Longueur des lectures (L)
 cat supplemental_data/Raw_stats_fastq.txt
 
#file                          format  type   num_seqs        sum_len  min_len  avg_len  max_len
#FASTQ/SRR10390685_1.fastq.gz  FASTQ   DNA   7,066,055  1,056,334,498       35    149.5      151
#FASTQ/SRR10390685_2.fastq.gz  FASTQ   DNA   7,066,055  1,062,807,718      130    150.4      151

## profondeur = N*L/G 
echo "(7066055*149.5 + 7066055*150.4)/4215606" | bc
#502

## Ecriture de la version de fastqc sur le fichier profondeur_seq_vs_genomeref.txt dans supplemental_data
echo "(7066055*149.5 + 7066055*150.4)/4215606" | bc > supplemental_data/profondeur_seq_vs_genomeref.txt

La profondeur de séquençage est de : 502 X.

Nettoyage des reads

Vous voulez maintenant nettoyer un peu vos lectures. Choisissez les paramètres de fastp [3] qui vous semblent adéquats et justifiez-les.

## Répertoire de travail
pwd
#/shared/projects/dubii2021/agodmer/Evaluation_M4_M5/DATA_analysis

## Chargement du module fastp
module avail fastp
module load fastp/0.20.0
## Ecriture de la version de fastp sur le fichier Version_tools.txt dans supplemental_data
#fastp --version >> ../supplemental_data/Version_tools.txt # ne fonctionne pas !
## Ecriture avec nano de la version du fichier

## Ressources sur le cluster
salloc --cpus-per-task=8 --mem=5G

## Lancement de la commande fastp avec les paramètres choisis
srun fastp --in1 FASTQ/SRR10390685_1.fastq.gz --in2 FASTQ/SRR10390685_2.fastq.gz --out1 CLEANING/SRR10390685_1_clean.fastq.gz --
out2 CLEANING/SRR10390685_2_clean.fastq.gz --html ./fastp.html --cut_mean_quality 30 --cut_window_size 8 --length_required 100 --cut_tail --json CLEANING/fastp.json

## Statistiques sur les read apres nettoyage avec fastp et écriture dans le fichier After_cleaning_stats_fastq.txt dans supplemental_data
srun seqkit stats --threads 1 CLEANING/*.fastq.gz > ../supplemental_data/After_cleaning_stats_fastq.txt

## Nombre de reads conservés
cat ../supplemental_data/After_cleaning_stats_fastq.txt

#file                                   format  type   num_seqs      sum_len  min_len  avg_len  max_len
#CLEANING/SRR10390685_1_clean.fastq.gz  FASTQ   DNA   6,777,048  996,891,051      100    147.1      151
#CLEANING/SRR10390685_2_clean.fastq.gz  FASTQ   DNA   6,777,048  990,442,597      100    146.1      151

## Calcul du nombre de reads conservés et ecriture dans le fichier Pctage_reads_non_conserves.txt (supplemental_data)
echo "scale=3;((7066055-6777048)/7006055)*100" | bc > ../supplemental_data/Pctage_reads_non_conserves.txt

Les paramètres suivants ont été choisis :

Parametre Valeur Explication
cut_mean_quality 30 Nettoyage des reads en fin de séquençage avec une qualité un peu plus basse
cut_window_size 8 Taille de la fênetre glissante pour le calcul de qualité moyenne des bases
length_required 100 Supprimer les reads ayant une taille inférieure à 100pb
cut_tail Enlever les adapateurs (non nécessaire, déjà réalisé)

Ces paramètres ont permis de conserver 6777048 reads pairés, soit une perte de 4.1% des reads bruts.

Alignement des reads sur le génome de référence

Maintenant, vous allez aligner ces reads nettoyés sur le génome de référence à l’aide de bwa [4] et samtools [5].

## Répertoire de travail
pwd
#/shared/projects/dubii2021/agodmer/Evaluation_M4_M5/DATA_analysis

## Chargement du module bwa
module avail bwa
module load bwa/0.7.17

## Il n'y a pas de commande "Version", écriture de la version du module dans le fichier Version_tool.txt avec l'éditeur de texte nano

## Création d'un FASTA index sur le génome de référence
srun bwa index REFERENCE_GENOME/GCF_000009045.1_ASM904v1_genomic.fna.gz

## Ressources du cluster
salloc --cpus-per-task=32 --mem=5G

## Mapping en utilisant la commande bwa mem
srun bwa mem REFERENCE_GENOME/GCF_000009045.1_ASM904v1_genomic.fna.gz CLEANING/SRR10390685_1_clean.fastq.gz CLEANING/SRR10390685_2_clean.fastq.gz -t 8 > MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1_genomic_reference.1.sam

## Chargement du module samtools
module avail samtools
module load samtools/1.10

## Ecriture de la version de samtools sur le fichier Version_tools.txt dans supplemental_dat
samtools --version >> ../supplemental_data/Version_tools.txt

## Ressources du cluster
salloc --cpus-per-task=8 --mem=5G

## Conversion du fichier SAM en BAM
samtools view MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1_genomic_reference.1.sam -b > MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1_genomic_reference.1.bam

## Trier le fichier BAM avec la commande samtools sort
srun samtools sort MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1_genomic_reference.1.sam -o MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1_genomic_reference.1.sort.bam

## Indexage du fichier BAM trié
srun samtools index MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1_genomic_reference.1.sort.bam

## Supression des fichiers inutiles
rm -f MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1_genomic_reference.1.sam SRR10390685_on_GCF_000009045.1_ASM904v1_genomic_reference.1.bam

Combien de reads ne sont pas mappés ?

## Répertoire de travail
pwd
#/shared/projects/dubii2021/agodmer/Evaluation_M4_M5/DATA_analysis

## Statistiques de MAPPING
srun samtools idxstats MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1_genomic_reference.1.sort.bam > ../supplemental_data/SRR10390685.sort.bam.idxstats
srun samtools flagstat MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1_genomic_reference.1.sort.bam > ../supplemental_data/SRR10390685.sort.bam.flagstat

## Calcul du nombre de reads non mappés

cat ../supplemental_data/SRR10390685.sort.bam.flagstat
#13571369 + 0 in total (QC-passed reads + QC-failed reads)
#0 + 0 secondary
#17273 + 0 supplementary
#0 + 0 duplicates
#12826829 + 0 mapped (94.51% : N/A)
#13554096 + 0 paired in sequencing
#6777048 + 0 read1
#6777048 + 0 read2
#12746420 + 0 properly paired (94.04% : N/A)
#12769290 + 0 with itself and mate mapped
#40266 + 0 singletons (0.30% : N/A)
#0 + 0 with mate mapped to a different chr
#0 + 0 with mate mapped to a different chr (mapQ>=5)

## Nombre de reads non mappés (total des reads : "mapped" + "singletons" + unmapped) ; unmapped = "total" - "mapped" - "singletons"
echo "13571369-12826829-40266" | bc
#704274

704274 reads ne sont pas mappés.

Croisement de données

Calculez le nombre de reads qui chevauchent avec au moins 50% de leur longueur le gène trmNF grâce à l’outil bedtools [6]:


## Répertoire de travail
pwd
#/shared/projects/dubii2021/agodmer/Evaluation_M4_M5/DATA_analysis

## Chargement du module bedtools
module avail bedtools
module load bedtools/2.29.2
bedtools --version >> ../supplemental_data/Version_tools.txt

## Extraction du gène d'intérêt à partir du fichier gff
srun zgrep trmNF REFERENCE_GENOME/GCF_000009045.1_ASM904v1_genomic.gff.gz | awk '$3=="gene"' > REFERENCE_GENOME/trmNF_gene.gff3

## Ressources du cluster
salloc --cpus-per-task=8 --mem=10G

## Croisemment des données avec l'option -f (nombre de reads qui chevauchent avec au moins 50% de leur longueur le gène trmNF)
srun bedtools intersect -a MAPPING/SRR10390685_on_GCF_000009045.1_ASM904v1_genomic_reference.1.sort.bam -b REFERENCE_GENOME/trmNF_gene.gff3 -sorted -f 0.5 > MAPPING/SRR10390685_on_trmNF_gene.bam

## Tri et indexage des reads chevauchants
srun samtools sort MAPPING/SRR10390685_on_trmNF_gene.bam -o MAPPING/SRR10390685_on_trmNF_gene.sort.bam
srun samtools index MAPPING/SRR10390685_on_trmNF_gene.sort.bam

## Nombre de read chevauchant le gène d'intérêt
srun samtools idxstats MAPPING/SRR10390685_on_trmNF_gene.sort.bam > ../supplemental_data/ SRR10390685_on_trmNF_gene_rapport.sort.bam.idxstats


## Visualisation du fichier
 cat ../supplemental_data/SRR10390685_on_trmNF_gene_rapport.sort.bam.idxstats

2801 reads chevauchent le gène d’intérêt.

Visualisation

Utilisez IGV [7] sous sa version en ligne pour visualiser les alignements sur le gène. Faites une capture d’écran du gène entier.

## Répertoire de travail
pwd
#/shared/projects/dubii2021/agodmer/Evaluation_M4_M5/DATA_analysis/REFERENCE_GENOME

## Dézipper le fichier fasta du génome de référence (peu volumineux)
gunzip GCF_000009045.1_ASM904v1_genomic.fna.gz

## Création d'un fichier FASTA index pour le génome de référence
samtools faidx GCF_000009045.1_ASM904v1_genomic.fna

## Téléchargement des fichiers avecla commande scp : 
#GCF_000009045.1_ASM904v1_genomic.fna
#GCF_000009045.1_ASM904v1_genomic.fna.fai
#SRR10390685_on_GCF_000009045.1_ASM904v1_genomic_reference.1.sort.bam
#SRR10390685_on_GCF_000009045.1_ASM904v1_genomic_reference.1.sort.bam.bai
#SRR10390685_on_trmNF_gene.sort.bam
#SRR10390685_on_trmNF_gene.sort.bam.bai

Capture d’écran du gène entier :
Capture d’écran du gène entier

References

1. toolkit NS. NCBI sra toolkit. NCBI, GitHub repository. 2019.

2. Andrews S. FastQC a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

3. Zhou Y, Chen Y, Chen S, Gu J. Fastp: An ultra-fast all-in-one fastq preprocessor. Bioinformatics. 2018;34:i884–90. doi:10.1093/bioinformatics/bty560.

4. Li H. Aligning sequence reads, clone sequences and assembly contigs with bwa-mem. arXiv preprint arXiv:13033997. 2013.

5. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and samtools. Bioinformatics. 2009;25:2078–9.

6. Quinlan AR, Hall IM. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.

7. Thorvaldsdóttir H, Robinson JT, Mesirov JP. Integrative genomics viewer (igv): High-performance genomics data visualization and exploration. Briefings in bioinformatics. 2013;14:178–92.

 

A work by Migale Bioinformatics Facility

https://migale.inrae.fr

Our two affiliations to cite us:

Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France

Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France